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Abstract 

In these lectures we give an overview of nonequilibrium stochastic systems. In par- 
ticular we discuss in detail two models, the asymmetric exclusion process and a 
ballistic reaction model, that illustrate many general features of nonequilibrium 
dynamics: for example coarsening dynamics and nonequilibrium phase transitions. 
As a secondary theme we shall show how a common mathematical structure, the 
(^-deformed harmonic oscillator algebra, serves to furnish exact results for both sys- 
tems. Thus the lectures also serve as a gentle introduction to things g-deformed. 
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1 Introduction 



In these lectures we explore the subject of nonequilibrium dynamics. Before 
getting into any kind of detail let us first establish what we mean by a nonequi- 
librium system. This is best done by taking stock of our understanding of an 
equilibrium system. Consider the Canonical (Boltzmann) distribution for a 
systems with configurations labelled C each with energy E{C): 

p(c) = ?5^tEEm (1) 
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where P — 1/ kT. The task is to calculate the partition function 



Z = ^exp(-/?E(C)), 



(2) 



c 



from which all thermodynamic properties, in principle, can be computed. The 
distribution (1) applies to systems in thermal equilibrium i.e. free to exchange 
energy with an environment at temperature T . It can easily be generalised to 
systems free to exchange particles, volume etc but it always relies on the 
concept of the system being at equilibrium with its environment. 

If one were interested in dynamics, for example to simulate the model on a 
computer, one might choose transition rates between configurations to satisfy 



where W{C' — > C) is the transition rate from configuration C to C. Condi- 
tion (3) is known as detailed balance and guarantees (under the assumption 
of ergodicity to be discussed below) that starting from some nonequilibrium 
initial condition the system will eventually reach the steady state of thermal 
equilibrium given by (1). We will discuss further this dynamical relaxation 
process and the properties of the steady state endowed by the detailed bal- 
ance condition in Section 3. For the moment we note that a system relaxing 
to thermal equilibrium is one realisation of a nonequilibrium system. In recent 
years such relaxation dynamics have been of special interest, for example, in 
the study of glassy dynamics whereby, on timescales realisable in experiment 
(or simulation), the system never reaches the equilibrium state and it is a very 
slowly evolving nonequilibrium state that is observed. This is sometimes re- 
ferred to as 'off- equilibrium' dynamics. Also let us mention the field of domain 
growth whereby an initially disordered state is quenched (reduced to a tem- 
perature below the critical temperature for the ordered phase) and relaxes to 
an ordered state through a process of coarsening of domains. The interesting 
physics lies in the scaling regime of the coarsening process which is observed 
before the equilibrium (ordered) state is reached. 

The other meaning of nonequilibrium refers to a system that reaches a steady 
state, but not a steady state of thermal equilibrium. Examples of such nonequi- 
librium steady states are given by driven systems with open boundaries where 
a mass current is driven through the system. Thus the system is driven by its 
environment rather being in thermal equilibrium with its environment. 

A pragmatic definition of a nonequilibrium system that encompasses all of 
the scenarios above is as a model defined by its dynamics rather than any 
energy function i. e. the configurations of the model are sampled through a local 
stochastic dynamics which a priori does not have to obey detailed balance. 



W{C' ^ C)e-^^(^') = W{C ^ C')e-'^-^(^) 



(3) 
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1.1 Structure of these notes 



These notes broadly follow the four lectures given at the summer school. In 
addition a tutorial class was held to explore points left as exercises in the 
lectures. In the present notes these exercises arc included in a self-contained 
form that should allow the reader to work through them without getting stuck 
or else leave them for another time and continue with the main text. The notes 
are structured as follows: in section 2 we give an overview of two simple models 
that we are mainly concerned with in these lectures. In section 3 we then set 
out the general theory of the type of stochastic model wc are interested in and 
point out the technical difficulties in calculating dynamical or even steady- 
state properties. Section 4 is an interlude in which we introduce, in a self- 
contained way, a mathematical tool — the gr-deformed harmonic oscillator — 
that will prove itself of use in the final two sections. In Section 5 we present 
the solution of the partially asymmetric exclusion process and amongst other 
things how the phase diagram (Figure 3) is generalised. In Section 6 we discuss 
the exact solution of a stochastic ballistic annihilation and coalescence model. 



2 Two simple models 

In this work we will focus on two exemplars of nonequilibrium systems: the 
partially asymmetric exclusion process and a particle reaction model. These 
models have been well studied over the years and a large body of knowledge 
has been built up [1]. We introduce the models at this point but will come 
back to these models in more detail in Sections 5 and 6 in which we summarise 
some recent analytical progress. 

2.1 Asymmetric exclusion process 
2.1.1 Model definition 

The asymmetric simple exclusion process (ASEP) is a very simple driven lat- 
tice gas with a hard core exclusion interaction [2]. Consider M particles on 
a one-dimensional lattice of length say. At each site of the lattice there is 
either one particle or an empty site (to be referred to as a vacancy or hole) — 
there is no multiple occupancy. 

The dynamics are defined as follows: during each time interval AT each parti- 
cle has probability AT of attempting a jump to its right and probability qAT 
of attempting a jump to its left; a jump can only succeed if the target site is 
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Fig. 1. Dynamics of the partially asymmetric exclusion process with open bound- 
aries. 

empty. 

In this work will be concerned with the case AT = dt ^ so that we obtain 
continuous time dynamics with particles attempting hops to the right with rate 
1 and hops to the left with rate q. In this limit no two particles will jump at 
the same time. The other limit of AT = 1 would correspond to fully parallel 
dynamics where all particles attempt hops at the same time. This type of 
dynamics is employed (for g = 0) in traffic flow modelling (see Schadschneider 
this volume). 

To complete the speciflcation of the dynamics we have to flx the boundary 
conditions. It turns out these are of great signiflcance. The following types of 
boundary conditions have been considered: 

Periodic In this case we identify site N + 1 with site 1. The particles then 
hop around a ring and the number of particles is conserved. 

Open boundaries In this case particles attempt to hop into site 1 with rate 
a. A particle at site leaves the lattice with rate f3. Thus the number of 
particles is not conserved at the boundaries. One can also understand these 
boundary conditions in terms of a site being a reservoir of particles with 
flxed density a and site + 1 being a reservoir with fixed density 1 — (3. 
We shall primarily be concerned with these boundary conditions because, 
as first pointed out by Krug [3], they can induce phase transitions. The 
dynamics is illustrated in Figure 1. 

Closed boundaries In this case the boundaries act as reflecting walls and 
particles can not enter or leave the lattice. Thus one has a zero current of 
particles in the steady state and the system obeys detailed balance. We shall 
not be interested in this case. 

Infinite system Finally we could consider our lattice of size A^ as a flnite 
segment of an inflnite system. 



2.1.2 Density and current in the ASEP 

The macroscopic properties of the ASEP in which we are interested are the 
steady-state current (flux of particles across the lattice) and density proflle 

(average occupancy of a lattice site). These can be expressed in terms of the 
binary variables {tj}, where = 1 if site i is occupied by a particle and 
Tj = if site i is empty, and which together completely specify a microscopic 
conflguration of the system. Then, the density at site i is deflned as {Ti{t)) 
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where the angular brackets denote an average over all histories of the stochastic 
dynamics. One can think of this average as being an average over an ensemble 
of systems all starting from the same initial configuration at time 0. 

Let us consider for the moment the totally asymmetric dynamics g = (i.e. no 
backward particle hops are permitted). One can use the 'indicator' variables 
{tj} intuitively to write down an equation for the evolution of the density. For 
example, if we choose a non-boundary site (so that we can avoid prescribing 
any particular boundary conditions) we obtain 

%^ = (r.-i(l-rO)-(r,(l-r,+0)- (4) 

Note that (rj_i(l — Tj)) gives the probability that site i — 1 is occupied and 
site i is empty. Thus, since particles hop forward with rate 1, first term on the 
right hand side gives the rate at which particles enter site i and the second 
term gives the rate at which they leave site i. 

If one needs more convincing of the argument used to obtain (4) consider what 
happens at site i in an infinitesimal interval dt: 



Ti{t + dt)^ < 



Ti{t) with probability 1 — 2dt 

Ti{t) + [l-Ti{t)\Ti^i{t) with probability dt ■ (5) 

ri(t)ri+i(t) with probability dt 



The first equation comes from the fact that with probability 1 — 2dt (dropping 
terms of 0(dt^)), neither of the sites i — 1 or i is updated and therefore Tj re- 
mains unchanged. The second and the third equations correspond to updating 
sites i — 1 and i respectively. If one averages (5) over the events which may 
occur between t and i -|- dt and all histories up to time t one obtains (4) . 

The same kind of reasoning allows one to write down straight away an equation 
for the evolution of (rjTj+i). 

^^^^^'^^^ = (Ti-l(l - Ti)Ti+i) - (TiTi+i(l - Ti+2)) (6) 

or for any other correlation function (rjTj- ■ ■ ■ ) . These equations are exact and 
give in principle the time evolution of any correlation function. However, the 
evolution (4) of (tj) requires the knowledge of (rjrj+i) which itself (6) requires 
the knowledge of (rj_irj+i) and (rj_irjTj+i) so that the problem is intrinsically 
an A^-body problem in the sense that the calculation of any correlation func- 
tion requires the knowledge of all the others. This is a situation quite common 
in equilibrium statistical mechanics where, although one can write relation- 
ships between different correlation functions, there is an infinite hierarchy of 
equations which in general makes the problem intractable. 
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For the case of open boundary conditions the number of particles in the system 
is not conserved. The evohition equations (4-6) arc then vahd everywhere in 
the bulk but they are modified at the boundary sites. For example (4) becomes 



d(ri) 



a((l-rO)-(ri(l-T2)) 



(7) 



dt 



dM 
dt 



{tn-i{1 - tn)) - P{tn) ■ 
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In the steady state, where the time derivatives of correlation functions are 
zero, (4,7,8) can be rewritten as a conserved current J 



J = a{{l - n)) = (ti(1 - r2)) = • • • (t,_i(1 - t,)) = • • • {tn-i{1 - tj,)) = /3(r^) . 

(9) 



These equations simply express that for the density to be stationary one must 
have the current into any site equal to the current out of that site. For example, 
a{{l — Ti)) is the probability that the leftmost site is empty multiplied by the 
rate at which particles are inserted onto it, and hence yields the current of 
particles entering the system. Similarly (3{tn) is the rate at which particles 
leave the system. 

2.1.3 Applications of the ASEP 

To motivate the study of this model let us consider a few applications. First 
we may consider the ASEP as a very simplistic model of traffic flow with 
no overtaking (more realistic generalisations are discussed by Schadschneider, 
this volume). Joining many open boundary systems together into a network 
forms a very simple model for city traffic. Such a network appears to represent 
faithfully real traffic phenomena such as jamming [4]. 

A major motivation for the study of the ASEP is its connection to interface 
dynamics and hence to the KPZ equation (noisy Burgers' equation). The KPZ 
equation is a stochastic non-linear partial differential equation that is central 
to much of modern statistical physics. In these lectures we do not have time 
to develop the theory of the KPZ equation, rather we refer the reader to [5]. 
Here we just make clear the connection to a particular interface growth model. 

The asymmetric exclusion model may be mapped exactly onto a model of a 

growing interface in (1 + 1) dimensions [6] as shown in Figure 2. The mapping 
is obtained by associating to each configuration {rj} of the particles a config- 
uration of an interface: a particle at a site corresponds to a downwards step 
of the interface height of one unit whereas a hole corresponds to an upward 



6 



step of one unit. The heights of the interface are thus defined by 

hi+i -hi = l-2Ti. (10) 

The dynamics of the asymmetric exclusion process in which a particle at site 
i — 1 may interchange position with a neighbouring hole at site i, corresponds 
to a growth at a site i which is a minimum of the interface height i.e. if 
hi{t) = hi_i{t) — 1 = — 1 then 

hi hiit) + 2 with rate 1 . (11) 

A growth event turns a minimum of the surface height into a maximum thus 
the ASEP maps onto what is known as a single-step growth model [6] meaning 
that the difference in heights of two neighbouring positions on the interface 
is always of magnitude one unit. Since whenever a particle hops forward the 
interface height increases by two units, the velocity v at which the interface 
grows is related to J the current in the asymmetric exclusion process by 

v^2J . (12) 



Periodic boundary conditions for the particle problem with M particles and 
N — M holes correspond to an interface satisfying /ij+Ar = hi + N — 2M, i.e. 
to helical boundary conditions with an average slope 1 — 2M/N. The case of 
open boundary conditions corresponds to special growth rules at the bound- 
aries. Because of this equivalence, several results obtained for the asymmetric 
exclusion process can be translated into exactly computable properties of the 
growing interface. 

One can also go on to list further applications [7]. Indeed early applications 
concerned biophysical problems such as single-filing constraint in transport 
across membranes [8] and the kinetics of biopolymerisation [9] . 



Fig. 2. The single-step growth model and its mapping onto the asymmetric exclusion 
process with q = 0. Note that the addition of a new particle to the surface (shown 
Ught-shaded) corresponds to the rightward hop of a particle in the ASEP. 
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(iii) 



(i) 



a 



Fig. 3. Phase diagram for the totally asymmetric exclusion process (g = 0). 
2.1.4 Phase diagram for q = 

Let us now discuss some interesting results that have been derived over the 
last decade for the open boundary system [10,11]. For the moment we present 
them without proof, although we shall recover some of them in the analysis of 
later sections. In Figure 3 we present the phase diagram (corresponding to the 
limit N ^ 00) ior the totally asymmetric exclusion process i.e. when q — 
as predicted within a mean- field approximation [9,12]. It should be noted that 
we consider the steady state in the limit 00, thus it is implicit that we 

have taken the limit t ^ 00 first. Note there are three phases — by phase it is 
meant a region in the phase diagram where the current has the same analytic 
form; phase boundaries divide regions where the current J and bulk density 
p (defined as the mean occupancy of a site near the centre of the lattice as 
N ^ 00) have different forms: 



(i) High density phase For /? < | and a> (3, 

J = /3(l-/3) 
(a) Low density phase For a < | and (3 > a, 

J — a{l — a) p — a 
(iii) Maximal current phase For a > | and /5 > |, 



T 1 1 



(13) 



(14) 



(15) 



Consider the high density phase. In this phase the exit rate is low and the bulk 
density is controlled by the right boundary. One can think of this as queue of 
traffic formed at a traffic hght that doesn't let many cars through at a time. 
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The low density phase is the opposite scenario where the left boundary controls 
the density through a low input rate of particles. One can think of this as cars 
being let onto an empty road only a few at a time. 

Finally in the maximal current phase the current of particles is saturated i. e. 
increasing a or f3 any further does not increase the current or change the bulk 
density. In this phase the density of particles decays from a its value at the 
left boundary to 1/2 its value in the bulk as x^^^'^ where x is the distance from 
the left boundary. 

This phase diagram of the model has many interesting features. Firstly we 
have both continuous (from high density or low density to maximal current) 
and discontinuous (from low density to high density) phase transitions even 
though it is a one-dimensional system. In the former transitions there are 
discontinuities in the second derivative of the current. In the latter although 
the current is continuous across the phase transition its first derivative and the 
bulk density jump. Along the transition line a = (3 < 1/2 one has coexistence 
between a region of the high density phase adjacent to the right boundary 
and a region of the low density phase adjacent to the left boundary. These 
results are very suggestive that the current J can be thought of as some 
free energy. We shall see later through an exact solution how this can be 
quantified. Also note that throughout the maximal current phase one has 
long-range, power-law correlations functions i.e. long-range correlations are 
generic. This contrasts with usual behaviour in equilibrium systems where 
power-law correlations are seen at non-generic, critical points (an exception 
is the Kostcrlitz-Thouless phase that may emerge at low temperature in two- 
dimensional equilibrium systems). 

2.2 Reacting particle systems and applications 

We now have a look at reaction-diffusion systems wherein particles can react 
with each other to form some by-product or more simply annihilate each other. 
Well-studied reaction systems are 

^ + ^^0 A + A^A. (16) 

The first reaction signifies that when two A particles meet they react and 
annihilate each other, the second that the two reagents coalesce. These systems 
model molecular reactions and the dynamics of laser induced quasiparticles 
known as excitons [13]. However the particles in this type of system can also 
represent composite objects such as aggregating traffic jams or, as we shall 
see, domain walls in phase-ordering kinetics. 

Consider first the A+A — > reaction. Clearly the steady state of such a system 
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is the uninteresting dead state where ah particles (except perhaps one) have 
been annihilated. The true interest lies in the late time regime wherein the 
system may enter a scaling regime. By this it is meant that the system is 
statistically invariant under rescaling by a typical length scale that depends 
on time. In this case the length scale can be taken as the typical distance 
between particles. 

For diffusing A particles, known results are that for d = 1 the typical dis- 
tance between particles grows as t^^^ or equivalently the density of a particles 
decreases as p ~ We shall recover this result in the next section. The 

result differs from a naive mean-field description obtained by assuming 

9p 2 

which would predict p ~ t~^. A similar density decay law holds for the 
A + A ^ A process and it has been shown that the two processes are ba- 
sically equivalent [14-17]. In addition to having a range of applications, dif- 
fusive reaction systems have served as prototypes for the development of a 
variety of theoretical tools including field-theoretic techniques [18] and the 
renormalisation group [19] and exact methods in one dimension [1,7]. 

One can also consider ballistic particle motion. In this class of models, parti- 
cles move deterministically with constant velocity and on meeting have some 
probability of coalescing into one particle or annihilating. A seminal paper 
by Elskens and Frisch [20] that considered the one- dimensional deterministic 
case where right moving and left moving particles always annihilate on con- 
tact, showed the decay depends on the initial densities of the two species. In 
particular, p ~ only if there are equal initial densities of the right-moving 
and left-moving particles. Ballistic models apply to chemical reactions when 
the inter-reactant distance is less than what would be the mean free path of 
particles. Such models may also be applied to various domain growth problems 
[21] and as we shall describe in section 6, the smoothing of a growing interface 
[22]. 



3 Theory of Stochastic Processes 



Having now given a rough overview of some of the interest in low dimensional 
nonequilibrium systems, the rest of these lectures will be aimed at a given 

a deeper understanding of how these phenomena come about in some of the 
simple models we have discussed. In this section we give a brief introduction 
to stochastic processes focussing on particular aspects that will be relevant to 
the models to be studied in detail in later sections. 
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We consider a system that can be in a finite number of microscopic config- 
urations and whose configuration changes according to transition rates. The 
simplest example is of a particle performing a random walk in continuous time 
on a one-dimensional lattice: the particle hops to the right with rate p and 
to the left with rate q. (The meaning of a 'rate' W is that in an infinitesimal 
time interval dt an event happens with probability W dt.) Since probability 
is conserved we can write a continuity equation for the probability which is 
referred to as 'the master equation' 

dP(x.t) 

' = pP{x - 1, t) + qP{x + 1, t) - (p + q)P{x, t) (17) 

The first two terms on the right hand side represent 'rates in' from other 
configurations {i.e. different positions of the particle); the last term represents 
the rate out of the given configuration. Note that the master equation is linear 
in the probabilities. 

For the case of the random walker, the master equation (17) can be solved 
analytically: it is a diffusion equation in discrete space but continuous time. 
If we consider a periodic lattice of A'^ sites (site A^-|- 1 is identified with site 1) 
the preferred method is to calculate the Fourier modes. Defining 

1 ^ 

P{k,t)^-J24PM with Zk^e'''''/'' (18) 

x=l 



yields from (17) 

dP{k,t) 
dt 

Thus 



XkP{k, t) where = p{zk - 1) + qiz^^ - 1) . (19) 



P{x,t)^^e^^'z^^P{k,0) with P(A;,0) = -^z^P(x,0). (20) 

k=l x=l 

The steady state corresponds to the mode k = N {\ = 0). The eigenvalues 
yield decay times Tk — 1/ |ReAfe| of each mode. For k finite and N large 

Afe ^ (p - q)^ - (p + q)^^ (21) 

Thus the equilibration time goes like r ~ where the dynamic exponent 
z = 2 as expected of a diffusion process. 



Exercise 1: Random walker on the Id lattice 

(a) Review the solution of equation (17) using Fourier modes described above. 
Show that if the initial condition is -P(x, 0) = bx^x^, the probability distribution at 
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time t is given by 



N 



P{x, t\xo, 0) = ^ XI ( ^'^^ ~ 



k=l 



27rik{x — xq] 
N 



where \k is given by (19) 

If you are unfamiliar with the discrete Fourier transform, you should first show that 
for l<k,e<N 



1 ^ k 



= Ok 



x=l 



and hence 



N N 



x=l 



k=l 



where Zk is as given by equation (18). 

(b) Consider now the same random walk problem but on an infinite rather than 
periodic lattice. Show from equation (17) that the generating function F{z,t) = 
EZ-ocPi^MQ/Pr^^^"" obeys 



dF{z,t) 

di 



[-{p + q) + Vpq{z + z-^)] F{z, t) 



and hence 



F{z, t) = F{z, 0) exp[-(p + q)t] exp[VPQ{z + z-^)t] 



(El.l) 



Given that the generating function of modified Bessel functions of the first kind 
In{x) is defined as 



X Inix)z"- = exp 

n=—oo 



{z + z~^)x 



show by expanding (El.l) that the general solution for the random walk problem 



is 



P{x,t) = eM-iP + Q)t] [ -) P{Wx-e{^VPQt) . (E1.2) 



Some of the properties of the modified Bessel functions /„(a;) can be gleaned by 
considering the special case of a symmetric random walker that begins at the origin, 
i.e. the case p = q = ^ and P{x,0) = dxfl. Then, P{x,t) = exp(— t)Ia;(i). Show 
then that the property /n(0) = 5n,o immediately follows, as does the fact that 
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I-n{t) = In{t)- Also, show from (17) that the Bessel functions satisfy the differential 
equation 



More generally we consider 'many-body' problems. For example a system of 
many particles on a lattice obeying stochastic dynamical rules. For concrete- 
ness we consider the simple example of the one-dimensional Ising model in 
which at each site i there is a spin Si — ±1 and periodic boundary conditions 
are imposed. 

As mentioned in the introduction the dynamics of an equilibrium model is 
usually chosen to satisfy detailed balance. In the case of the Ising model under 
spin-flip dynamics, whereby each spin can flip with a certain rate according 
to the directions of the neighbouring spins, the detailed balance condition on 
the spin-flip rates of spin Si at site i reads 

w[-S^^Si) = + S,^,)S.] (22) 

(We have taken E — —J2iSiSi+i). In Glauber dynamics (22) is satisfied by 
choosing 

W{S, ^ -S,) = ^-S^t.nmS,_, + S,^,)] ^23) 

as is easily verified. Then, in the T = limit the spin fiips for a down-spin 
{Si = — 1) at site i occur with the rates indicated in Figure 4 and similarly for 
up-spins at T = 0. Note the final transition where the spin flips against both 
its neighbours is forbidden. Thus domains of aligned spins can not break up. 
The domain walls between aligned domains i.e. neighbouring pairs of spins 1 1 
or 1 1 make random walks on the bonds of the original lattice and annihilate 
when they meet. As noted above, this process of diffusing 'particles' (here the 
domain walls) that annihilate on meeting is sometimes denoted A-l-A 0. The 
solution for the mean domain length has been known for some time [23,14] — 
basically it can be solved by considering the random walk performed by the 
length of a domain (distance between neighbouring particles) . 

To sec the solution for the mean domain size within the kinetic Ising model 
[24] consider the master equation which can be written 

^^^^ = -J2PiC,t)W{S, ^ -S,) + Y,P{Q,t)W{-S, ^ S,) (24) 
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Fig. 4. T = dynamics of the Ising model under Glauber dynamics. Also indicated 
is the interpretation of domain walls as annihilating particles. 

where C is a configuration of the Ising spins and Cj is the same configuration 
with spin flipped. The number of domain walls is given by ^(1 — SiSi+i) /2 

i 

SO to study the coarsening process, in which domain walls are eliminated 
and aligned domains increase in size, one considers the equal time correlation 

function Ck{t) — {j;^'^Si{t)Si+k{t)) where the average indicates an average 

i 

over initial conditions and possible histories of the stochastic dynamics. It can 
be shown that obeys a discrete diffusion equation 

^^^^^^ ^ -2C{k,t) + C{k + l,t) + C{k-l,t) for k>0 (25) 
with boundary condition C{0,t) = 1. 

The fact that we obtain a closed set of equations for two point correlation 
functions is a very happy property that allows one to solve easily for these 
correlation functions — basically (25) is an equation of the same form as (17) 
but with a boundary condition corresponding to a source of walkers at site 
0: see [25,26] for details. Note that the method generahses to finite tempera- 
ture but not to higher dimensions. Actually higher point correlation functions 
are more difficult to calculate; nevertheless it has been shown that the one- 
dimensional Glauber model can be turned into a free fermion problem which 
then implies, in principle at least, that any correlation function can be calcu- 
lated [27]. 

Remarkably, the full domain size probability distribution has been explicitly 
calculated for the one-dimensional g-state Potts model, which includes the 
Ising model as the case q = 2 [28]. (Note that a calculation such as Exercise 2 
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only calculates the first moment of the domain sizes). The dynamics of domain 
walls in the g-state Potts model, represented as A particles, corresponds to 
diffusion and, on meeting, reaction according to 



A + A 



A with probability (g — 2)/(g — 1) 

(26) 

with probability l/{q — 1) 



To understand this equation note that when a domain is climated each of 
the neighbouring domains may be in one of g — 1 states so that they have 
probability — 1) of being in the same state and therefore coalescing. Thus 
for g = 2 we recover A + A — > and for g — > oo (infinite state Potts model) 
we recover A + A ^ A. 



Exercise 2: Time dependence of the Id Ising model 

(a) The mean magnetisation M{t) in the Id Ising model is defined as 

AT N 

^(*) = N ^^^^^ ^ TV ^ ^ ^^-(^^ 

J=l {C} j=l 

where N is the number of spins on the lattice, {C} is the set of all 2^ possible spin 
combinations and Si{C) is the value of the spin at site i, in configuration C. Show 
for general single spin-flip dynamics as described by (24) that the magnetisation 
obeys the differential equation 

{C} i=i 



Under Glauber dynamics defined by equation (23) show that 

^ = -(l-tanh2/3)M(t) 

and hence that M{t) = M(0) exp(— [1 — tanh2(3]t). Also comment on the finite- and 
zero-temperature behaviour (recalling that /3 is inverse temperature). 

Hint: you will need to argue that the relation tanh ^{Si-i+Si+i) = ^ - Si tanh 2/3 
holds for any (periodic) combination of spins. 

(b) Show that the mean density of domain walls noit) in the Id Ising model can be 
expressed using the spin correlation function C{k, t) as no (t) = l[l-C{l,t)]. Under 
zero-temperature Glauber dynamics, the evolution of C{k,t) is described by the set 
of difference equations (25) which has the stationary solution C{k) = IV A; > 0. 
Verify this solution, and explain its physical meaning. 
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To obtain the time-dependence of C{k,t) we introduce a function e{k,t) through 
C{k,t) = 1 — e{k, i). The time-evolution is governed by equation (17) with P{k^ t) = 
e{k,t), p = q = 1 and subject to the boundary condition e(0,t) = 0\/t. 

This boundary condition can be satisfied through an appropriate choice of the initial 
conditions e{k, 0) in the 'unphysical' region k < (this is equivalent to using the 
Method of Images [29]). Using the general solution (El. 2) and the properties of the 
modified Bessel functions derived in exercise 1, show that e(0, t) = Vt if e{—k, 0) = 
e(A;,0). Hence show that 

oo 

C{k,t) = 1 - exp(-2t) ^ [Ik-ei2t) - Ik+i{2t)] 

e=i 

is the solution for the spin correlation function that has C(0,t) = IVf and with 
the initial spin orientations uncorrelated, i.e. C{k,Q) = Skfi. Show that this implies 
that for these initial conditions, noit) = ^ exp(— 2t)[7o(2t) -|- Ii{2t)]. 

Finally, using the asymptotic form of the modified Bessel functions 



exp(x) 

71^ 



TTX 



^ k\(8x) 

k=l ^ ' j=l 

show that the long-time behaviour of noit) is noit) (47rt)~^/^ -|- 0(i~^/^). 



In order to formalise stochastic processes in general we consider a master 
equation 

^E!^A^ p{c')w{c' ^c)- Yl p{c)w{c^c') (27) 

where W{C ^ C) is the rate of transition from configuration C to C . It is 
useful to combine the transition rates into a single matrix 



M{C,C') = W{C' ^C) for C'^C (28) 
M{C,C) = -Y,W{C ^C) (29) 



c 

so that we can write 

dP(C,t) 



= ^M(C,C')P(C',i). (30) 



Now let us introduce a 'bra-kef notation. Let the vector of probabilities at 
time t (one component for each configuration) be \P)t- Thus the probability 
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of the system being in configuration C at time t is (C|-P)t and we can write 
the master equation as 



^ = M\P), . (31) 

Thus the master equation looks something hke a Schrodinger equation with 
'imaginary time' and B. replaced by —M (this is why M is sometimes referred 
to as a Hamiltonian). However it should be recalled that in quantum mechanics 
it is the modulus squared of the components \ P)t that are probabilities whereas 
here the components of \P)t are probabilities. See [30,7,31], for example, for 
more details. 

The formal solution of the master equation (31) is 

|P)t = exp(MO|P)o (32) 

where |P)o gives the initial conditions. However, a more useful way of dealing 
with the problem is to consider the eigenvectors of the transition matrix M. 
The first thing to note is that since M is not Hermitian, in general, it has 
different left and right eigenvectors 

M|0a) = A|0a) (33) 
(V'a|M = A(V'a| . (34) 

These satisfy 'bi-orthogonality' 

(V'a|0;.) = V> (35) 
and in principle (assuming M can be diagonalised) we may write 

1 = EI'^a)(V'a| (36) 

A 

M = Y^\\<p^){^^\ . (37) 

A 

Then we have using (32), (36) and (33) 

= E —T- = E^^p(^^)I'^a)(^a|P)o . (38) 



Note that the real parts of all eigenvalues A are negative except for one eigen- 
value which is zero. The reason for this is often referred to as the Perron- 
Frobenius theorem [32]. This theorem holds for a matrix with non- negative 
entries where a sufficiently large power of the matrix has all positive entries. It 
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states that the maximum eigenvalue is non-degenerate and the corresponding 
left and right eigenvectors have positive entries. 

In our context one can relate the spectrum of M to that of non-negative matrix 
M obtained by adding some suitable multiple (say minus the diagonal element 
with the largest magnitude) of the identity matrix to M. The assumption 
that a sufficiently large power of M has all positive entries is equivalent to 
the assumption that the (finite) system is ergodic i.e. any configuration can 
be reached from any other after sufficient time. In that case the system will 
have a unique steady state so that the zero eigenvalue cannot be degenerate 
and the steady state eigenvector has positive elements which correspond to 
probabilities. 

In principle, all probabilities at all times can be found if one can diagonalise 
the matrix M. Unfortunately this is only possible in a few isolated cases. 
Generally this corresponds to many-particle systems where the particles do 
not interact or certain 'integrable' systems that can be solved via the Bethe 
ansatz. 

The Bethe ansatz was originally used in the context of a spin-chain model for 
magnetism [33] and is suitable for a one- dimensional lattice system in which 
the number of particles is conserved. It is basically a generalisation of the 
Fourier transform. Let us briefly discuss it for the case of the ASEP. The idea 
is to guess the following form for eigenvectors 

0(xi, . . . , xm) = E A(Q)zt'''''z',''''' ■ ■ . (39) 
Q 

in which Xi is the position of particle i on a lattice, and 1 < aji < 2:2 < ^3 . . . < 
Xm < N with the lattice size once again. The sum is over all permutations Q 
of (1,2,... , M) and A{Q) is some amplitude. One uses the master equation to 
determine relations that must be satisfled by (f>{xi, . . . ,xm) in its unphysical 
regions, e.g. where two particles occupy the same site and Xi — Xj+i for some i. 
After some manipulation, one arrives at a set of nonlinear equations that must 
be simultaneously satisfied by all the 'momenta' Zi. An early account of this 
method in the context of stochastic processes is given in [34]. The technique 
is also described in [7,35] . However even for the so-called 'integrable' systems 
where the Bethe ansatz works it is actually very difficult to extract the explicit 
eigenvectors. 

Of course, it may not be possible to solve a model exactly, and in this situation 
one often turns to approximation schemes. For example, one can try a 'short 
time' expansion — we simply truncate the power series defining the exponential 
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of M in (32) after m terms 



-E^I^>o- (40) 



n=0 ^■ 

The finite number of terms m on the right hand side can be calculated. From 
the series one can try to estimate, for example, the dynamical exponent [36]. 

A project less ambitious than solving the full dynamics is to calculate the 
steady state probabilities (the eigenvector with eigenvalue 0). We denote the 
steady state by dropping the subscript t so that \P) is the eigenvector with 
eigenvalue zero 

^ = M|P) = , (41) 

and its components P{C) are the steady state probabilities. Thus one has the 
balance conditions 

J2 W{C C')P{C) = ^(^' ^ C)P{C) . (42) 
c c 



Let us review some simple ways to satisfy (42). First of all one has the case of 
detailed balance where 

W{C C')P{C) = W{C' C)P{C') . (43) 

However to check this condition i. e. to check whether detailed balance is satis- 
fied we first need to know the steady state! However there is an equivalent set 
of conditions called 'Kolmogorov criteria' which state that detailed balance is 
obeyed if and only if the transition rates satisfy 



W{Ci ^ C2)W{C2 ^ C3) . . . W{Cn ^ Ci) 

= W{Ci ^ Cn)W{Cn ^ Cn-l) ■ ■ ■ W {C2 ^ Ci) (44) 

for every possible cycle in configuration space Ci,C2, . . . ,C„,Ci. 

That (43) implies (44) is very simple to show. To show that the converse is 
also true is most simply achieved [37,38] by using (43) to determine P{C). 
Then, one finds that P{C) is unique if and only if (44) is satisfied. 

Let us consider some of the implications of the detailed balance condition. 

First note that when one has detailed balance then in the steady state there 
is no net flow of probability between any two configurations. Since there is 
no fiow of probability there is nothing to distinguish the forwards direction in 
time from backwards direction. Therefore running the systems backwards in 
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time will not change, for example, any two time correlation functions. Thus 
the system is said to be reversible. 

Now consider the simple problem of a single random walker on a periodic 
one-dimensional lattice, studied at the beginning of this section. The steady 
state is given by equal probability of the walker occupying each site. However, 
one does not have detailed balance as defined above unless the hop rates are 
symmetric p = q. But one knows that the only effect of p 7^ g is to impose 
a drift. Clearly with a drift the system is not reversible. But the system is 
reversible under simultaneous reversal of time and parity {i.e. direction): the 
parity operation results in defining the image of position * — N + 1 — X. 

Then trivially 

W{x x+l)P{x) = W{N-x N-x + l)P{N-x) (45) 

since in the steady state P{x) is constant in space and both transition rates 
are equal to p. 

More generally an extended detailed balance condition can pertain as defined 
in [39]: 

W{C C')P{C) = W{C'* C*)P{C'*) (46) 

where C* is the image configuration of C under a reversible mapping. One 
requires (C*)* —C and 

c c 
Again one can construct an equivalent condition as follows [40] 

W(Ci ^ C2)W{C2 ^ C3) . . . W{Cn ^ Ci) 

= w{ct ^ c:)w{c: ^ c:_,) . . . w{c; ^ c*) (48) 

for every possible cycle in configuration space Ci, C2, C3, . . . ,Cn along with the 
additional condition (47). 

The condition (46) is known in the probabilistic literature as dynamic re- 
versibility [38]. More recently it has appeared in the physics literature under 
the name pairwise balance [41]. It holds for example in the ASEP on a ring 
where the image configuration is obtained by the parity operation. 
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3. 1 Formal solution for the steady state 



In this subsection our aim is to derive the general solution for the steady 
state (41) for an arbitrary transition rate matrix M. Although the solution is 
too general to be of any practical use it is instructive to see how a quantity 
analogous to a partition function naturally arises. 

As it stands, the set of linear equations (41) is underdetermined because the 
determinant of M is zero (one of its eigenvalues is zero). However, we know that 
a probability distribution must be properly normalised, so consider instead the 
equation 



(49) 



in which M*^*) is the matrix obtained from —M (the minus sign is introduced 
for convenience below) by replacing the i*^ row with all ones, i.e. the vector 
(l,l,l,...l). The ket vector |i) is the basis vector corresponding to configu- 
ration i. This set of equations can now be solved and we do so using Cramer's 
rule. This states that the probability of the system being in configuration Cj 
in the steady state is 



det M(^;^) 
det M(^) 



(50) 



in which the matrix M^^'^^ is obtained from Af*^') by replacing column j with 
a vector with a 1 at position j and zeros everywhere else. For simplicity we 
choose i = j, and to be clear we write out M^'^^ explicitly. 



1 



-M. 



111 



-M„,_i -M, 



1 



-M, 



(51) 



We note that the determinant det M^^'^^ is equal to the determinant of the 
matrix obtained from M = —M by removing row j and column j. This deter- 
minant is called a cofactor of M and will be written fj = f{Cj). Furthermore, 
we observe that the determinant in the denominator of (50) is the sum of all 
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these cofactors. That is 

n 

Z = T.fj (52) 

so that 

P{Cj) = f . (53) 

We thus identify the cofactor fj with the steady-state weight of configuration Cj 
and Z as a normalisation for the stochastic process. In an cquiUbrium system, 
Z is related to the partition function Z(.q = J2c^^Pi~(^^i^))- However, the 
normahsation Z obtained using the above procedure is not guaranteed to be 
identical to Zeq due to the possibility of common factors present in all the 
weights fj and hence Z. With this understood, (52) extends the notion of a 
partition function to nonequilibrium systems. 



3.2 Relationship of the steady-state normalisation to eigenvalues of the transition- 
rate matrix 



We now note an interesting and, as far as we are aware, little-known result that 
relates the normalisation Z (52) to the eigenvalues of the transition matrix. 
The relation reads 

z = n i-^j) (54) 

in which the product is over the eigenvalues of the transition rate M except 
the zero eigenvalue. 

This result is obtained by expanding the characteristic polynomial det(AlL— M) 
in powers of A. One finds [42] 

det(Al - M) = (-)" I det M + (-)"A J2 fj + 0(^^) | (55) 

with fj as defined above. Given that det M — and Z — J2j fj we find that 

^ det(Al-M) n,(A-A,) -n- , . ^ 

Z = lim ^ = lim ^ = TT -A, 56 

A-.0 A A 

where the last step follows because one, and only one, eigenvalue of M is zero 
(at least if the process is ergodic). 
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3.3 Nonequilibrium phase transitions 



We now consider the existence of phase transitions in a stochastic process. As 
in equihbrium systems, we identify phase transitions through discontinuities 
in macroscopic quantities or their derivatives. Furthermore, at a first-order 
(discontinuous) transition we expect phase coexistence and at a continuous 
transition divergences in length and time scales. 

In terms of eigenvalues of the transition matrix M, a phase transition would 
be indicated by the clustering of one or more eigenvalues towards zero. For 
example as a discontinuous transition is approached one would expect long- 
lived 'excited' states since the decay time of an eigenstate of M is inversely 
proportional to the real part of the corresponding eigenvalue. These long- 
hved states can be thought of as metastable states and the approach of an 
eigenvalue to zero in the thermodynamic (infinite system-size) limit suggests 
phase coexistence. Near a continuous transition one has diverging time scales 
(as manifested by e.g. a power-law decay of an autocorrelation function) so 
one expects a continuous spectrum of eigenvalues whose real parts are close 
to zero. 

Bearing this in mind, we note that equation (54) implies that the steady-state 
normalisation Z also approaches zero at a phase transition in the thermody- 
namic limit. Such a scenario is analogous to Yang-Lee theory of equilibrium 
phase transitions, in which it has been rigorously shown that zeros of the grand 
canonical partition function in the complex-fugacity plane converge towards a 
point on the positive real axis that indicates the value of the fugacity at which 
a phase transition occurs. Similar properties are also displayed by the zeros of 
the canonical partition function in the complex-temperature plane (these are 
also known as Fisher zeros). Hence, it would seem hkely from (54) and the 
above considerations that zeros of a steady-state nonequilibrium normalisation 
in the complex plane of some control parameter would also converge to the 
real axis at a nonequilibrium phase transition. This would also imply nonan- 
alyticities in a macroscopic variable averaged with respect to the steady-state 
probability distribution function. 

This idea has been pursued in the context of an asymmetric exclusion process 
[43]. In this case, the steady-state normalisation was solved for its zeros in the 
complex-fugacity plane and indeed they were found to approach the real axis 
for values of the model parameters that corresponded to a phase transition. 
Additionally, a similar approach has been tried in the context of directed 
percolation [44], a further lattice model that exhibits a nonequilibrium phase 
transition (see e.g. [31] for a review). 

This concludes our overview of stochastic processes. We will shortly move on 
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to the specific models we are going to study. As mentioned in the introduction 
the 'g-deformed Harmonic Oscillator' plays a central role in the solution of 
both these models. 



4 Intermezzo: What is a g'-deformed Hcirmonic Oscillator? 



First let us recall the quantum harmonic oscillator well known from under- 
graduate physics. The Hamiltonian is 

H-\a' + \f (57) 

where we have lazily set all constants m, % tt etc to unity and where " denotes 
an operator. One way to treat the problem is to introduce two new operators 
a, through 

1 i 

^"Tf^^^^^-* "Tl^^"^^-* ■ ^^^^ 

Then inserting these definitions into the commutation relation [x,p] = 
yields 

ad^ - a^a = 1 . (59) 
Also the Hamiltonian (57) becomes 

^ = 0^0+^. (60) 

The operators d, are, of course, raising and lowering operators. More gen- 
erally they are bosonic operators, the defining property of which is the com- 
mutation relation (59). In the number (or energy) basis they obey 

a|n) = 1^ ~ 1) d^\n) = Vn -|- 1 |n -|- 1) for n>0. (61) 

Thus in this basis d^d\n) = n |n) and the Hamiltonian is diagonal. The number 
n labels the energy eigenstatcs. The quantised energy excitations obey Bose 
statistics and the operators (61) create and annihilate bosons respectively. 

The Hermite polynomials come into play when we consider the projection of 
an energy eigenstate onto the position eigenbasis defined by 

x\x) — x\x) . (62) 

Then 

= e-^'i/„(x) (63) 
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where Hn{x) is the Hermite polynomial of degree n. 

We now consider a 'deformation' of (59) by introducing a parameter q [45] 

aa) - qa)a = 1 . (64) 

Observe that taking q = 1 recovers bosonic operators whereas q = —1 gives 
fermionic operators. Thus it was originally hoped that the gr-deformed opera- 
tors might describe some new excitations interpolating between fermions and 
bosons. 

The operators a and a) now operate on basis vectors \n) (with n = 0,l,2,...) 
as follows: 



V 1 - ? 



— 1^ ~ 1) ■ (66) 

Thus a)d, which is the analogue of the energy operator, is diagonal in this 
basis 

1 — g" 

a^a\n) — |n) . (67) 

The eigenstates of the 'position operator' x — a-\-a} are given by 

x\x) = \x) . (68) 

\/l — q 

Finally, we identify the gf-Hermite polynomials as 

H^{x\q) = {x\n) . (69) 



4-1 Properties of q-Hermite polynomials 

Using (65,66) and (69) it is straightforward to derive a recurrence relation 



^Jl - q-+^H^+i{x\q) - 2xH^{x\q) + ^1 - q-H^_M<l) = . (70) 
This can be compared to the usual recurrence relation for Hermite polynomials 
i/„+i(x) - 2xHn{x) + 2nH„-i{x) = . (71) 
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To make the connection between g'-Hermite and ordinary Hermite polynomials 
we would like the recurrence relations to coincide in the limit q ^ 1. The 
prescription is 

Hn{x\q) - ((l -g)/2r^' ff^(^(i_g)/2 ^) (72) 

Q)n 

in which both the Hermite polynomial Hn{x\q) and the independent variable 
X are transformed. The factors of \/2 that appear can be traced back to the 
fact that we set x = (a + a^) / \pl when doing the quantum mechanics, but we 
take £ = a + aMn the g-deformed case. 

Explicit formulae for can be found using a generating function technique, 
the details of which differ slightly depending on whether g < 1 or g > 1. Here 
we discuss in detail the case q<\. 

First we define a generating function G{x^ A) for the g-Hermite polynomials 

_ oo \ n 

G(x,A) = E^==W- (73) 

n=o \J{q]q)n 

where have introduced a 'g-factorial' notation defined through 



n— 1 

(a;g)„= n(l-0 (74) 

(a;g)o = l. (75) 

We later encounter products of these factorials for which we use a standard 
shorthand [46]: 

(a, 6, c, . . . ; g)„ = (a; g)„(6; g)„(c; g)„ . . . . (76) 

The g-factorial in (73) is introduced for convenience as will now become ap- 
parent. 

We obtain a functional relation for (5 (a;. A) by multiplying both sides of equa- 



tion (70) by A"/ \J (g; g)„ and performing the required summations: 

G{x, qX) = (A^ - 2Xx + l)G{x, A) . (77) 

It is convenient to parameterise the 'position' eigenstates by an angle 9 where 
X — cos 9 and < ^ < 27r, and to replace G{x, A) by G{9, A) (a function of 9). 
Then (77) becomes 
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Since q < 1 we can iterate to obtain 

G{e,x) 



(79) 



where we have used the g- factorial notation (76). Note that ^(6*, 0) is just 
(^|0) which we are free to set to 1. 

The infinite product l/{x;q)ao has a well-known and easy to verify series 
representation [46] valid for \x\ < 1, q < 1 



1 



oo 



(x;q). 



E 



(t, q), 



(80) 



from which, with a little effort, we may extract the form of {0\n). Expanding 
both sides of (79) in A and comparing coefficients we find 



{0\n) 



E 



i{n-2k)0 



where 



is the g-deformed binomial defined as 

{Q; Q)n 



(81) 



(q;q)n-k(q;q)k 



(82) 



when < A; < n and zero otherwise. In the limit q 1 the g-binomial 
coefficient is equal to the conventional version familiar from combinatorics. 

Important properties of the g-Hermite polynomials arc their orthogonality and 
completeness. It can be shown [46] that the set of g-Hermite polynomials are 
orthogonal with respect to a weight function i'{9). That is 



{n\9)u{e){e\m) = Sn, 

J 

where the weight function is given by 

... (g,e^-^e'^-^;g)oo 
^^^^ ^ 2^ 



(83) 



(84) 



Completeness implies we can form a representation of the identity matrix: 

r<ie\e)v{9)(e\^ \ . (85) 

Jo 



It is less well known that, in a similar manner to the above exposition, one can 
obtain explicit forms for g-Hermite polynomials when q> 1 [47]. It turns out 
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that they can be obtained from (81) by making the substitution 9 — > 7r/2 — iu 
i.e. one replaces cos 6 by isinhu and the parameter u runs from —oo to oo. 
The weight function that orthogonahses these polynomials is now 

1 1 

''^''^ ~ h^ (g-l,-g-le2«,-g-le-2-g-l)^ " ^^^^ 



Exercise 3: Coherent states, the exponential function and their g- 
defor mat ions 

For the 'traditional' quantum mechanical raising and lowering operators that satisfy 
the commutation relation (59) construct the vector |/x) = exp(/ia^)|0). Show that 
this is a coherent state, i. e. an eigenvector of the lowering operator a with eigenvalue 
II. 

Now consider a g-coherent state |/i) satisfying the eigenvalue equation {l — q)^/'^a\fi) = 
li\fi) with a and the g-deformcd raising and lowering operators that obey (64). 
By solving this eigenvalue equation, show that the g-coherent state |//) may be 
expressed as 

l/x) =exp,([l-g]i/W)|0) 
where the g-deformation of the exponential function is defined through 

exp,(x) = — . 

n=0 

Show that this scries converges always for q > 1, for < 1 when q < 1 and that 
limg^i expg([l - q]x) = exp(a;). 

Verify also that when |g| < 1 and expg{x) converges, it can also be written as 

^"^^^"^ = (^ • 

Hint: this is most easily achieved by checking that both the series and product 
representations of expg(x) satisfy the functional relation expq(gx) = (1 — x) exp^(a;) 
and that they are numerically equal for a special value of a; = 0. 

Finally, show that 

oo 

{e\i,) = J2mn){n\fi)=G{e,fi) 

n=0 

where G{9,iJ,) is the generating function of g-Hermite polynomials (73,79). 
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5 Steady State Solution of Partially Asymmetric Exclusion Process 



We now consider in more detail the partially asymmetric exclusion process 
(PASEP) introduced in Section 2.1. Recall that we discussed some of the 
properties of the special case where the reverse hop rate q was set to zero 
and in particular we presented a phase diagram for the model (Figure 3). We 
now describe how one obtains such a phase diagram exactly for general q. 
We do this by relating the steady-state probability distribution of the model 
to products of matrices which can be mapped onto the gr-deformed harmonic 
oscillator ladder operators. 



5. 1 The matrix product formulation and its quadratic algebra 

In this section we review the matrix approach to finding the steady state of the 
model. The approach has been the subject of lectures at a previous summer 
school in this series [35] and the reader is referred there for further details. 
Here the bare essentials of the method are outlined. 

Consider first a configuration of particles C and its steady-state probabil- 
ity P{C). We use as an ansatz for P{C) an ordered product of matrices 
X1X2 . . . Ajv where Xi = D if site i is occupied and Xi = E if it is empty. To 
obtain a probability (a scalar value) from this matrix product, we employ two 
vectors {W\ and \V) in the following way: 



Note that here the matrices, bra and ket vectors are in an auxiliary space 
and are not related to the space of configurations and the probability ket 
vector we considered earlier. The factor Z^v is included to ensure that P{C) is 
properly normalised. This latter quantity, analogous to a partition function as 
was discussed in Section 3, has the following simple matrix expression through 
which a new matrix C = D+E is defined: 



Note that if D and E do not commute P{C) is a function of both the number 

and position of particles on the lattice, as expected for a non-trivial steady 
state. The algebraic properties of the matrices can be deduced from the master 
equation for the process [10]. It can be shown that sufficient conditions for 
equation (87) to hold are 



P{C) 




(87) 



Zn = {W\{D + Ef\V) = {WIC^'IV) . 
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DE - qED = D + E (89) 
a{W\E={W\ (90) 
PD\V) = \V). (91) 

In this formulation steady state correlation functions are easily expressed. For 
example, the mean occupation number (density) of site i may be written as 

{W\C'-'DC^-'\V) 
in} = y ■ (92) 

To get a feel for why the conditions (89-91) give the correct steady state 
consider the current J of particles between sites i and i + 1. 

J = {nil - Ti+i)) - q{{l - ri)n+,) (93) 

In the matrix product formulation the current becomes 

J ^ {W\C'-\DE - qED)C''-'-'\V) ^^^^ 

Now using (89) and the fact that C — D + E we find 

j^iWjC^^Z^^ (95) 

We see that, as is required in the steady state, the current is independent of 
the bond chosen. 

To see that the current into site one reduces to the same expression we use 
(90) 



a{W\EC^-'\V) Z^_, 
J = q;(1-Ti) = (96) 

Also the current out of site N reduces to the same expression using (91) 

J = P{tn) = ^ = (97) 

Z/7V 

Note that the fact that all the above expressions for the current are equivalent 
is a necessary condition for the matrix formulation to be correct. A sufficient 
condition is to check that the master equation (41) is satisfied for all 2^ con- 
figurations. The algebra (89-91) allows one to do this systematically [10,48]. 

Our task now is to evaluate the matrix products in the above expressions for 
Zjv, J and Tj by applying the rules (89-91). 
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In [10] the case q = was treated by using (89) repeatedly to 'normal-order' 
matrix products: that is, to obtain an equivalent sum of products in which all 
E matrices appear to the left of any D matrices. For example, consider powers 
of C 



C" = l 

C^ = D + E 

C^^D^ + E^ + {1 + q)ED + D + E 



N N-n 



n=0 m=0 

Then finding a scalar value would be straightforward using (90) and (91): 

Zj, = {WlC'lV) = J2 aN,n,ma-''/3-''' . (98) 

n,m 

The difficulty with this approach lies in the combinatorial problem of finding 
the coefficients aN^n,m- This was solved for g = in [10]. However, the solution 
(and the problem of actually calculating the current and correlation functions) 
for arbitrary q remained open for some time although the phase diagram was 
conjectured [49]. Recently the problem has been overcome by making the 
connection with the g-Hermite polynomials discussed in Section 4 [50,51]. 

5.2 Connection with the q-deformed harmonic oscillator 

To make the connection with the g-deformed harmonic oscillator, let us define 



1 1 

+ , a (99) 

E^^—^ , ^ gt (100) 
1-g ^JY^q ^ ^ 

One finds that the algebraic rules (89-91) reduce to 



aa^-qa}a^\ (101) 

l^\a^ - -^^L=(H^| where w = - 1 (102) 

a\V) = ^ \V) where v = - 1 . (103) 

V 1 - ? P 
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Thus D and E are related to the g-bosonic operators discussed in Section 4; 
such a relationship was first pointed out in [49]. Also we see that {W\ and \ V) 
are eigenvectors of the g-bosonic operators. Such eigenvectors are known as 
coherent states (see exercise 3). In the oscillator's "energy" eigenbasis {|n)} 
(61) one can check that they have components 

{n\V) = ^== {W\n) = ^== . (104) 

For the moment we consider v < 1, w < 1 so that the vectors (104) are 
normalisable. 

We introduced earlier a matrix C which appears in the expressions for the 
mean particle density and current. We now see that this matrix can be written 
as a linear combination of the identity 1 and the "position" operator x — a+a^: 

C ^D + E^^—l+ } X . (105) 
1 - g VI - ? 

As we saw in Section 4 the eigenstates of the oscillator in the co-ordinate 
representation are continuous g-Hermite polynomials. Clearly, the eigenvectors 
of C are the same as those for x and therefore knowledge of them permits 
diagonalisation of C. 

Let us illustrate how to apply the procedure to obtain an expression for the 
normalisation Zjv when q<\. First we insert a complete set of states into the 
expression for the normahsation (88): 

Z^^ rAQv{Q)l^W\C^\Q)l,Q\V) . (106) 
h 

By design, the matrix C is acting on its eigenvectors, so using (105) and 

^ , ^. 2 cos ^ , , , 

x\d) = (107) 

we obtain 



\\qv{Q)I^W\B)['^-^^^^^-^\ {Q\V). (108) 
\ 1 — q J 



Now we know the form of {W\ and \V) in the {\n)} basis (104). Therefore we 
insert a complete set of the basis vectors in (108) to find 

oo oo n 

{9\V) = Y.m{n\V) = E ^==(^1^) • (109) 
n=o n=Q^{q;q)n 
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Observe that the final sum in this equation is nothing but (73), the generating 
function of the g-Hermite polynomials! Thus, when |v| < 1 and |w| < 1, we 
may write 

{e\V) = G{9, v) {W\9) = G{9, w) (110) 
where G{9,X) is given by expression (79). 

Putting all this together, we arrive at an exact integral form for the normali- 
sation 

^'d^i/(^)[2(l + cos^)]^G'(^,w)G'(^,^;) (111) 
which, written out more fully and using the notation (76), reads 

^^-^^l,r^J y^d^[2(i + cos^)] (^,.,^^,-.,^^,.,^^,-.,.^)^- 

(112) 

When \v\ > 1 or |w| > 1 equation (111) is not well-defined because G{9,\) 
does not converge when |A| > 1. Rather than finding a representation of the 
quadratic algebra (89-91) that does not suffer from this problem, one can 
simply analytically continue the integral (112) to obtain Z^v when \v\ or 1^1 
takes on a value greater than one. Since a contour integral is defined by the 
residues it contains, to analytically continue the integral one simply has to 
follow the poles of the integrand as they move in and out of the original 
integration contour (see exercise 4). 

One can also apply the procedure of the previous section to find an integral 
representation of the normalisation for the case ol q> 1. The only difference 
is that we must use the q> 1 polynomials. It turns out that this amounts to 
substituting cos^ by isinhw with the limits on u now running from — oo to 
-|-oo and using the weight function (86) 

(1 \^ oo 
I rduviu)\2il + isiTihu)fGiu,w)Giu,v). (113) 
1 — qJ J-oo 

It is possible to derive an alternative expression for which takes the form 
of a finite sum rather than an integral and is valid for all values of the model 
parameters. Such an expression implies the solution of the combinatorial prob- 
lem of reordering operators D,E under the rule (89). In [51] the formula (eq. 
55 in that paper) was derived by explicitly evaluating the integral (111) using 
further properties of g-Hermite polynomials and their generating functions. 
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5.3 Phase diagram 



The phase diagram for the partially asymmetric exclusion process is now ob- 
tained by calculating the current through the relation J = Zjv-i/^at and the 
above exact expressions for the normalisation Z^. 

For the case of forward bias {q < 1) in which particle hops are biased from 
the left boundary (where particles are inserted) to the right boundary (where 
they are removed) one expects a nonvanishing current in the thermodynamic 
limit N ^ oo. This is calculated by applying the saddle-point method to the 
integral (112) — see exercise 4. It turns out that for large N 

{W\C''\V) ~ a'^N-^ (114) 

Then in the thermodynamic limit J — > 1/a. Phase transitions arise from 
different asymptotic forms of the current due to the analytic continuation 
(see exercise 4) for the integral (112). Ultimately one finds the expressions 
presented in Table 1 and thence the phase diagram shown in Figure 5. 

Note that the current as given by (96) is a ratio of two partition functions of 
size N — 1 and N. In equilibrium statistical mechanics (in the large N hmit) 
this ratio is equivalent to the fugacity z (within the grand canonical ensemble). 
Recalling that z = and the chemical potential fi is equivalent to the Gibbs 
free energy per particle, gives credence to our claim in section 2.1.4 that the 
current J acts as a free energy for the system. 

We note that this phase diagram has a very similar form to that of Figure 3 
which applied for the special case of total asymmetry q = 0. That is, one finds 
the same three phases, namely a (i) high density, (ii) low density and (iii) 
maximal current phase. As discussed in Section 2.1 the transition between the 
low and high density phases (i) and (ii) is first order and those to the maximal 
current phase (iii) are second order. Note from Table 1 that, for example, the 



Phase 


Region 


Current J 


(i) 




P{l-q- (5) 
l-q 


(ii) 


oi < I3>a 


a{\ — q — a) 
l-q 


(iii) 




l-q 
4 



Table 1 

The N ^ oo forms of the particle current in the forward biased phases. 
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Fig. 5. The phase diagram for the current in the forward-bias regime of the PASEP. 

first derivative of the current is discontinuous at the transition between the 
low and high density phases. 

Similar methods can be used to find the density profiles in the PASEP: the 
details of the calculations and results are presented in [52] . It turns out that 
the low- and high- density phases subdivide into three subphases for nonzero 
q less than one [52]. In each of the low density subphases, the bulk density is 
p = «/(! — g); however what distinguishes the subphases is that the decay of 
the density profile from the right boundary takes a different form in each. The 
same is true of high-density phase, except that p — (1 — — gr) and that 

the density decay is at the left boundary. 



Exercise 4: Evaluation of the current in the PASEP with g < 1 

Show that the integral representation of the normalisation for < 1, equation (112), 
can be recast as a contour integral 

Zn = j> dzexp{Nh{z))g{z) (E4.1) 

where C is the unit circle centred on the origin, 

h{z) =ln{2 + z + z-^)-ln{l-q) and g{z) - ^ iq,z^z^-q)^ 



Airiz {vz,v/z,wz,w/z;q)c 



The large number of singularities in the integrand at the origin makes this a difficult 
integral to evaluate exactly using, e.g., the residue theorem. Instead, we shall use 
this contour integral to obtain the currents in the PASEP for large system size. 

First, we will apply the saddle-point method to the above integral which is valid for 
the range \v\ < l,\w\ < 1. The idea is to deform the contour of integration such that 
it passes through a saddle point in Re[/i(2:)] along the path of steepest descent (this 
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is also the line where lm[h(z)] = const). Show that this saddle point is at zq = 1 
and that the path of steepest descent (at least near zq) is parallel to the imaginary- 
axis. Also, show that g{zQ) = g'izo) = 0. 

Insert Taylor expansions of h{z) and g{z) (to second order) about the saddle point 
zo into (E4.1) to obtain 

ZN = eMNh{zo)) £dzew (^-^Nh" izo)iz - zof^ ^^{z - zof . 

You should now convince yourself that the deformation of the integration contour 
to pass through zq along the path of steepest descent in Re[/i(z)] is equivalent to 
making the substitution z = zq + ix with x running from — oo to +00. Hence show 
that the resulting Gaussian integral will give rise to an expression of the form 

Zn ~ Ki exp(iV/i(zo))iV-3/2 . 

Note: to save time, do not evaluate the constant Ki (although you may like to check 
that it is real and positive). 

Now show that the current in the PASEP for |i;| < 1, < 1 and large N behaves 
as 

Zn-1 1 l-q 

J = ^ = . 

Zpf exp h{zo) 4 



There now remains the question of what to do when one of the control parameters, 
say V, has a magnitude greater than 1. When \v\ < 1, the contour C encloses the 

poles of g(z) dX z = v, qv, q^v, . . . but excludes those at z = 1/v, q/v, q^ /v, Asv 

is increased, the pole at z = v moves outside C and that at 2; = 1/v moves inside. To 
obtain the analytic continuation of Zjv, one must add to the saddle-point expression 
for the integral over C above the residue at z = v and subtract that at z = 1/v. 
Show that this yields 

Zn ~ Ki exp(Ar/i(zo))iV-3/2 + K2 ex.p{Nh{v)) . 

Again try to avoid explicitly calculating the positive constant K2. Verify that h{v) > 
h{zo) and thus that, for sufficiently large N, Z^ ~ exp{Nh{v)). Hence show that 
for 1 < v < 1/q 

Zn-1 1 P{l-q-p) 

Zpf exp h{v) 1 — q 

Now convince yourself that as v is increased yet further, the contribution from the 
additional poles that have entered/exited the unit circle is smaller than that from 
those at z = v,l/v. Consider also the effect of increasing it; to a value less than v 
and show that when w > v,w > 1 

1 a{l — q — a) 

exp h{w) 1 — q 
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Lattice site — ^ 

Fig. 6. Proposed form of the density profile in the reverse-bias regime of the PASEP. 

Finally, use the information learned from this analysis to construct the phase dia- 
gram (Figure 5) for the current in the a-(5 plane. 



For g > 1 one has a reverse bias regime [51] where the boundary conditions 
force a current through the system against the direction of the reverse bias. 
This gives rise to a new phase in which the current decreases exponentially 
with the length of the system as J ~ g^^/^. This phase is of interest in 
the context of 'backbend dynamics' [53,54] where, for example, a fluid in a 
permeable medium has to traverse a pore oriented against the direction of 
gravity. 

One can predict the form of the density profile in the reverse bias regime as 
follows. We expect a high-density region at the left edge of the lattice, towards 
which the particles are biased but cannot escape, and a low-density region at 
the right edge of the lattice. Thus the density profile is a Fermi-like function as 
illustrated in Figure 6. In order that the profile is stationary the typical time 
for a particle to diffuse (against the bias) from the edge of the high-density 
region to the right boundary must be equal to the typical time of a hole to 
diffuse from the edge of the low-density region to the left boundary. Since 
the diffusion rates of both particles and holes are the same, the boundary 
between the high- and low- density regions must be in the centre of the lattice, 
as illustrated in Figure 6. 

5.4 Generalisations 

Here we give a brief overview of other related exclusion processes. Firstly the 
features of the phase diagram for q <1 (Figure 5) appear robust for stochastic 
one-dimensional driven systems and can be predicted from considerations of 
domain wall dynamics [55,56]. This has been confirmed through the exact 
solution of the open boundary system with fully parallel dynamics [57,58]. 
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However when the bulk dynamics become deterministic the maximal current 
phase is lost. 

For the case of symmetric exclusion q — 1 (no bulk drive) one does not have 
phase transitions; the current is easy to calculate and decays linearly with sys- 
tem size [59-61]. Intuitively one can understand this as diffusion between two 
particle reservoirs of different density. The diffusion current will be propor- 
tional to the concentration gradient which decreases as Thus although 
the bulk dynamics is symmetric the boundary conditions force a weak current 
of particle through the system. 

For g = 1 a recent work [62] has determined the probability of observing a 
particular coarse-grained density profile p{x). In an equilibrium system, this 
probability is given by exp(— A^jF[p]) in which jF[p] is a local functional that 
expresses the free energy difference between the equilibrium profile and p. 
From the matrix product approach a form for a functional that plays a similar 
role in the symmetric exclusion process was found in [62] and is shown to 
be nonlocal. For the particular boundary conditions where the dynamics sat- 
isfy detailed balance, J-'[p] reduces to the free energy functional known from 
equilibrium statistical physics and so the nonlocal functional would seem to 
generalise a well-known equilibrium concept. 

An interesting kind of boundary-induced phase transition, manifesting spon- 
taneous symmetry breaking, is found when the asymmetric exclusion process 
is generalised to two oppositely moving species of particle: one species is in- 
jected at the left, moves rightwards and exits at the right; the other species is 
injected at the right, moves leftwards and exits at the left [63]. Intuitively one 
can picture the system as a narrow road bridge: cars moving in opposite direc- 
tions can pass each other but cannot occupy the same space. The model has a 
left-right symmetry when the injection rates and exit rates for the two species 
of particles are symmetric. However for low exit rates this symmetry is 
broken and the lattice is dominated by one of the species at any given time. 
This implies that the short time averages of currents and bulk densities of the 
two species of particles are no longer equal. Over longer times the system flips 
between the two symmetry-related states. In the /3 ^ limit the mean flip 
time between the two states has been calculated analytically and shown to di- 
verge exponentially with system size [64]. Thus the 'bridge' model provided a 
first example of spontaneous symmetry breaking in a one dimensional system. 

In the models discussed so far the open boundaries can be thought of as 
inhomogeneities where the order parameter (particle density) is not conserved. 

Inhomogeneities which conserve the order parameter can be considered on a 
periodic system. Indeed a single defect bond on the lattice (through which 
particles hop more slowly) is sufficient to cause the system to separate into 
two macroscopic regions of different densities [65] : a high density region which 
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can be thought of as a traffic jam behind the defect and a low density region 
in front of the defect. Here the presence of the drive appears necessary for 
the defect to induce the phase separation. Moving defects {i.e. particles with 
dynamics different from that of the others) have also been considered and exact 
solutions obtained [66-69] . For the simple case of a slowly moving particle the 
phenomenon of a queue of particles forming behind it has been shown to be 
analogous to Bose condensation [70,71]. 



6 Stochastic Ballistic Annihilation and Coalescence 

In Section 2.2 we introduced two reaction systems, namely annihilation {A + 
yl — > 0) and coalescence {A + A^A).As noted in Section 2.2, it is well known 
that annihilation and coalescence reactions are equivalent when the reactant 
motion is diffusive [14,15] and hence lead to the same t~^^'^ density decay in 
one dimension. 

We now consider in detail the distinct case of ballistic reactant motion by 
which it is meant that particles move with some fixed velocity. We study a 
class of models that comprise an arbitrary combination of annihilation and co- 
alescence of particles with two (conserved) velocities and stochastic reaction 
dynamics and which was solved in [73]. As we shall see, an exact solution is 
possible by virtue of a matrix product method in which the matrices involved 
can be written in terms of the raising and lowering operators of the ^-deformed 
harmonic oscillator. A related, but distinct, ballistic reaction model (not con- 
tained within the class discussed here) had also revealed a connection between 
ballistic reaction systems and g-deformed operator algebras [74]. 



We now define the class of models to be considered. At time t — reac- 
tants are placed randomly on a line. Each particle is assigned a velocity +c 
(right-moving) or — c (left-moving) with probability and /l = 1 — /r respec- 
tively. Particles move ballistically until two collide, at which point one of four 
outcomes follows, see Figure 7: the particles pass through each other with 
probability q; the particles coalesce into a left (right) moving particle with 
probability prj^ {pVr}] the particles annihilate with probability p(l — rjL — rju). 
Here p — 1 — q is the probability that some reaction occurs. 

As an example of an application of this model consider the identification of 




Fig. 7. The possible reactions and their probabilities. 
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right- and left-moving particles with the edges of terraces as illustrated in 
Figure 8. If new particles are added to the system in such a way that they 
only stick to the sides, and the rate of particle addition is taken to infinity, one 
obtains ballistic motion of terrace edges. When two edges meet, a terrace is 
completed which corresponds to annihilation of approaching particles. Hence 
the relevant parameters of the general annihilation coalescence model are rjL — 
rjR = 0. The scattering reaction, which occurs with rate q, corresponds to the 
possibility of a new terrace being formed when two edges meet. 

The case of deterministic annihilation q — tjl — rjR — was studied by 
Elskens and Frisch [20]. Here particles always annihilate on contact. This case 
was found to exhibit a density decay that depends on time as t^^/^, but only 
if the initial densities of the two particle species (left- and right-moving) are 
the same. In the following analysis of the more general class of reaction sys- 
tems, we find that such a result persists if two 'effective' initial densities (to 
be defined below) are equal. Furthermore, we will see that the introduction 
of stochasticity into the reaction dynamics i.e. the parameter q, induces a 
transition in the density decay form. 

Our aim is to calculate the density decay. To do this we shall consider without 

loss of generality a right moving test particle as illustrated in figure 9. We 
define Pg (t) as the probability that the test particle survives up to a time t 
respectively. From figure 9 one can see that the initial spacing of the particles 
on the line does not affect the sequence of possible reactions for any given 
particle, in particular for the test particle (we return to this point later) . Also 
note that after a given time t, the test particle may only have interacted with 
the particles initially placed within a distance X = 2ct (and to the right) 
of the chosen particle. These two facts imply that the survival probability can 
be expressed in terms of two independent functions. The first is G{N; X), the 
probability that initially there were exactly N particles in a region of size 
X — 2ct. The second is F[N), the probability that the test particle survives 
reactions with the N particles initially to its right, and depends only on the 




Fig. 8. (i) Mapping of a surface growth model to a particle reaction system through 
the identification of upward and downward steps to left- and right-moving parti- 
cles. Addition of new particles that stick only to the sides of the surface cause the 
particles to move, (ii) The correspondence of particle annihilation and scattering to 
the elimination and nucleation of a terrace respectively. 
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Space 



Fig. 9. A configuration and set of trajectories and reactions for a test particle 
(shaded and bold line) encountering a string of iV = 10 particles. Note how changing 
the spacing between, for example, the fifth and sixth particles (indicated by dotted 
lines) alters the time sequence of the reactions but not the final survival probability. 

sequence of the particles. Explicitly, 

oo 

P^s''\x)=J2G{N;X)F{N). (115) 

N=0 

Thus the problem is reduced to two separate combinatorial problems of cal- 
culating G{N;X) and F{N). Note that the choice of G{N;X) allows one to 
consider a model defined on a lattice or on the real line. In the present work we 
assume particles are initially placed on a line with nearest-neighbour distances 
chosen independently from an exponential distribution with unit mean. Then 
we can use a well-known result [29] 

G{N;X) = X^e-^/N\ . (116) 

We now show that the second problem (calculation of F{N)) can be formulated 
within a matrix product approach. 

As an example consider a test particle encountering the string of reactants de- 
picted in Figure 9. We claim that the probability of the test particle surviving 
through this string may be written as 

{W\RRRLRLLLLR\V) (117) 

where R, L are matrices and {W\, \ V) vectors with scalar product (Vr|V") = 1. 
Thus we write, in order, a matrix R (L) for each right (left) moving particle 
in the initial string. 

The conditions for an expression such as (117) to hold for an arbitrary string 
are, in fact, rather intuitive 



41 



RL = qLR + p{r]LL + + [1 - 77L - m]) 
{W\L={W\{q + pr]R) 
R\V) = \V). 



(118) 
(119) 
(120) 



The condition (118) just echoes the reactions that occur in Figure 7 i.e. after 
an interaction between a right-moving and left-moving particle there are four 
possible outcomes (see Figure 7) corresponding to the four terms on the right 
hand side of (118) with probabihties given by the respective coefficients. Using 
(118), any initial matrix product such as (117) can be reduced to a sum of 
terms of the form {W\L^ R^\V) corresponding to all possible final states ensuing 
from the initial string and with coefficients equal to the probabilities of each 
final state. These final states give the possible sequences of particles that the 
right-moving test particle will encounter. The test particle will survive such 
a final state and pass through the s left-moving particles with probability 
[q + pri^y. The condition (120) ensures that this probability is obtained for 
each possible final state i.e. 



Finally the condition (119), along with = 1, ensures that once a right- 

moving particle passes through all the left-moving particles in the string it no 
longer plays a role. Thus (121) becomes 



Note that the reason for a matrix product formulation is different from the 
PASEP. There the steady state probability of any configuration, for arbitrary 
system size, could be written as a matrix product. Here the probability of 
a particle surviving some arbitrary sequence of particles can be written as a 
matrix product. 

The above approach relies on an important property of the system which is 
invariance of a reaction sequence with respect to changes of initial particle 
spacings. To understand this, consider again Figure 9. By altering the initial 
spacings of the particles, the absolute times at which trajectories intersect 
and reactions may occur (if the reactants have survived) may be altered. For 
example, by increasing the spacing between the fifth and sixth particles, the 
trajectories of the third and fourth particles can be made to intersect first. 
However as we have already seen, for any particle, the order of intersections 
it encounters does not change and so the final states and probabilities are 
invariant. This invariance is manifested in the matrix product by the fact that 
the order in which we use the reduction rule (118) is unimportant i.e. matrix 
multiplication is associative. Thus, it is the invariance with respect to initial 
spacings that allows the system to be solved by using a product of matrices. 



{W\VR'\V) = {q + pr}Ry{W\R'\V) . 



(121) 



{W\UR'\V) = {q+pr,Ry . 



(122) 
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We now proceed to evaluate F{N). Averaging (117) over all initial strings 
of length N, recalling that //? and fi are the probabilities that a particle is 
assigned velocity +c, — c respectively, yields 

F{N) = {W\{hL + fRRf\V) . (123) 

To evaluate F{N) we first write 



R^— a + 77L and a^ + r)R. (124) 

JR JL 

In these equations we have introduced two new parameters 

/fl = /i?(l-r/L) and /2 = /4l-r?^) (125) 
which we call effective initial particle densities and whose ratio 




(126) 



turns out to be an important parameter of the model. 

One can verify from (118-120) that a, satisfy a g-deformed harmonic oscil- 
lator algebra 



aa' — qaJa = 1 — q (127) 
i^W\a^ = {W\q/^ (128) 
a\V) = ^\V). (129) 

Thus (123) can be written as 

N 

\V) . (130) 

The vectors (W\ and \V) are eigenvectors of a^,a, and as in Section 5, they 
can be explicitly calculated. 

At this point one can see that we have precisely the same mathematical struc- 
ture as we did when solving for the steady state of the PASEP — compare 
(127-129) with (101-103). We just need to carry out the same procedure of 
diagonalising a matrix that is linear combination of the identity operator and 
the position operator of a g-deformed harmonic oscillator. Rather than repeat 
the details here we go directly to a discussion of the phase diagram. 



F{N) = {W\ 
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Exercise 5: Deterministic ballistic annihilation and random walks 



The Elskens-Frisch model [20] of deterministic ballistic annihilation is a special case 
of the stochastic model described above that has the parameters q = 'r]L = VR — 
0. For the case where particles are initially placed on all sites of a lattice with 
equal probability of right- and left-moving particles, write down all possible initial 
conditions that would allow a right-moving test particle to survive through one, 
two and three sites. Hence relate the survival probability of the test particle after N 
sites to the probability of a random walker not returning to the origin after + 1 
steps. Hint: devise a criterion for the survival of the test particle that involves the 
relative number of left- and right-moving particles for a given initial condition. 

Construct the matrices R and L and the vectors {W\ and \V) and explain how the 
resulting expression for F{N) (with fL = fR = 5) is related to the random walk 
problem. 



6.1 Phase diagram 

We consider the long-time density decay Q[t) which is obtained from the 
survival probabilities Pg^\t), Ps^\t) for right and left moving test particles 
through 



{Pg (t) can be deduced from P^ (t) through the left-right symmetry of the 
model.) 

The matrix product calculation for F{N) implies that the density decays as 



In this equation ^00 = ~ x) is the fraction of particles remaining once no 
more reactions are possible {i.e. all particles moving in the same direction). 
We now list the results the residual density Qoo when x < 1- (Results for x > 1 
one be deduced using the left- right symmetry of the model.) 

For x< 



Q{t) = fR- 



■pr{t)+fLPr{t). 



(131) 



Q{t) = ^00 + A(t) . 



(132) 




exp [-2 (1 - q) ifl - rjq) ct] 
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Fig. 10. Phase diagram for the model of stochastic balhstic annihilation and coa- 
lescence. The density decay form depends on the ratio of effective densities x 
the stochasticity parameter q. 
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R 



1 



exp 



-^ixffi-Mrct 



4 

oo 



4v27r(g^, q/y/x] q) 
( 1 



2 

oo 
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We now consider the phase diagram — figure 10. This is spanned by the val- 
ues of X (the ratio of effective densities x = fRlfl) ^"^^ ? (^^e stochasticity 
parameter): % parameterises all information concerning the asymmetry be- 
tween the right and left moving particles whereas q parameterises the level 
of stochasticity. Thus there is universality of ballistic annihilation and co- 
alescence: for a generic choice of rj^^ 7]^ defining a particular annihilation- 
coalescence model, the same four decay regimes are found by varying the 
initial densities or stochasticity parameter q. In this way the universality can 
be considered as a 'law of corresponding states'. 
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The line X = 1 (equal effective densities) is non-generic since a single, power 
law, decay regime is found. The decay does not depend on the stochasticity q. 
This special line was found by Elskens and Frisch [20] for the case of determin- 
istic annihilation {til = Vr = Q = ^) ^^'^ we have shown that this special line 
is present for all combinations of annihilation and coalescence. This phase can 
be understood through the picture of [20]. Density fluctuations in the initial 
conditions lead to trains of left- and right-moving particles: in a length ~ t 
the excess particle number is ~ t^^"^ which yields the density decay. At 
long times, the train size is large and so a particle in one train encounters 
many particles in the other and will eventually react making the parameter q 
irrelevant. 

The second phase found in [20] (labelled the Elskens-Frisch phase in the dia- 
gram) is seen to persist for nonzero q and in this phase the two particle species 
decay at equal rates. The two new phases that arise in the full model for x ^ 
{i.e. as a consequence of randomness in the reaction dynamics) have the con- 
trasting property that two particle species decay at unequal rates leaving a 
non-zero population of left-moving particles. A simple example of non-equal 
decays is the case rjR — 0,r]L — I {x = Q)- Then left-moving particles do 
not decay but simply absorb the right-moving particles with probability 1 — q 
giving g = /^e"^*^^"'^-''^^'^*. Our results show that, in general, increasing q leads 
to a non-trivial transition at x = to a regime where the two species have 
different decay forms. 



7 Conclusion 

In these lectures we have tried to give a flavour of the collective phenomena 
exhibited by simple low-dimensional systems with stochastic dynamics. Specif- 
ically we have seen that the partially asymmetric exclusion process exhibits a 
number of phase transitions (of both first and second order) and long-range 
correlations even in one dimension. We have also shown how similar mathe- 
matical techniques can be applied in the distinct context of a particle reaction 
system to obtain an exact solution. 

We hope that the reader with a background in equilibrium statistical physics 
is pleasantly surprised by the diversity of phenomena that emerge in low- 
dimensional nonequilibrium systems. We have also endeavoured to illustrate 
how novel analytical techniques are being developed for these systems. Through 
the inclusion of background material and the exercises, we hope also to have 
inspired confidence in the reader to approach the literature and make an active 
contribution to this fertile field. 

Acknowledgments: We would like to thank our collaborators F. Colaiori, F. 
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